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We introduce a novel variance-reducing Monte Carlo al- 
gorithm for accurate determination of autocorrelation times. 
We apply this method to two-dimensional Ising systems with 
sizes up to 15x 15, using single-spin flip dynamics, random site 
selection and transition probabilities according to the heat- 
bath method. From a finite-size scaling analysis of these au- 
tocorrelation times, the dynamic critical exponent z is deter- 
mined as 2; = 2.1665 (12). 



The onset of criticality is marked by a divergence of 
both the correlation length ^ and the correlation time r. 
While the former divergence yields singularities in static 
quantities, the latter manifests itself notably as critical 
slowing down. To describe dynamic scaling properties, 
only one exponent is required in addition to the static ex- 
ponents. This dynamic exponent z links the divergences 
of length and time scales: r ~ In our computation 
of z we exploit that for a finite system ^ is limited by 
the system size L, so that r ~ at the incipient critical 
point. 

In this Letter, we focus on the two-dimensional Ising 
model with Glauber-like dynamics. Values quoted in the 
literature for z vary vastly, from z = 1.7toz = 2.7[0, 
but recent computations seem to be converging towards 
the value reported here. Finally, results are beginning to 
emerge of precision sufficient for sensitive tests of funda- 
mental issues such as universality. 

The numerical method introduced in this Letter is re- 
lated to Monte Carlo methods used to compute eigenval- 
ues of Hamiltonians of discrete or continuous quantum 
systems ||,^ and transfer matrices of statistical mechan- 
ical systems Q]. In particular, the current method is suit- 
able to obtain more than one eigenvalue by adaptation 
of the diffusion Monte Carlo algorithm of Ref. ||] . 

To compute the autocorrelation time of small L x L 
lattices we exploit the following properties of the single- 
spin-flip Markov (or stochastic) matrix P Q. It oper- 
ates in the linear space of all spin configurations and its 
largest eigenvalue equals unity. The corresponding right 
eigenvector contains the Boltzmann weights of the spin 
configurations; the left eigenvector is constant, reflect- 
ing probability conservation. The correlation time tl (in 



units of one flip per spin, i.e., L'^ single-spin flips) is de- 
termined by the second-largest eigenvalue A^: 
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For a system symmetric under spin inversion, the corre- 
sponding eigenvector is expected to be antisymmetric. 

We used two methods to compute Xl- exact, nu- 
merical computation for L < 5, and Monte Carlo for 
4 < L < 15. The exact method used the conjugate gra- 
dient algorithm Q and the symmetries of periodic sys- 
tems. This calculation resembles that in Ref. [||, but 
currently we realize Glauber-like dynamics using heat- 
bath or Yang transition probabilities and random site 
selection. 

The Monte Carlo method used a stochastic form of 
the power method, as follows. A spin configuration s 
with energy E(s) has a probability 



exp(-S(s)/fcr) _ ^b(s)2 
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where Z is the partition function. The element P(s'\s) of 
the Markov matrix is the probability of a single-spin-flip 
transition from s to s'. Since P satisfies detailed balance. 
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is symmetric. For an arbitrary trial state |/) an effective 



eigenvalue A^*^ is defined by 
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where (•)/ is the expectation value in the state |/). In the 
limit t oo the effective eigenvalue converges generically 
to the dominant eigenvalue allowed by the symmetry of 
I/). The convergence is exponential in the time lag t. 

Given a trial state |/), standard Monte Carlo suffices 
to compute the right-hand side of Eq. (^), i.e., the de- 
nominator of Eq. (^ 

EE (/|P*|/) = 

J2 f{st+l)P{st + ll\st) ■ ■ ■ P(S2|S1)/(S1) = 
si,---,st+i 
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\'(/'B(si)V'B(st+l)/ p' 

is an autocorrelation; /(s) = and {■)p denotes the 

average with respect to the probabiUty 

Pist+i\st)---Pis2\si)Msi)VZ (6) 

of finding a configuration si in equihbrium and subse- 
quent transitions to configurations S2 through st+i. 
Similarly, the numerator of Eq. (^) 

jj(t) ^ (/|pt+i|/) = 
E f{st+i)P{st+i\st) ■ ■ ■ P(si|.so)/Gso) = 

2 \ VB(Si)tpB[St+l) / p 

is a cross-correlation, where the 'configurational eigen- 
value' Al(s) of spin configuration s is defined as 

Finally, with Eqs. (|) and (0), one has A^'^ = i/(*)/7V(*) 
for the effective eigenvalue. 

In practice, i/*^'^ and TV^*) are estimated by conven- 
tional Monte Carlo methods. As usual, these estimators 
involve time averaging of stochastic variables. Thus, on 
the right of Eqs. (||) and is replaced by St'+i-i 

{i — 1, . . . ,t), and the Monte Carlo average is taken over 
an appropriately chosen subset of times t' after thermal 
equilibration. 

In principle, one could choose / = mipB, where m is 
the magnetization. In that case, the above method re- 
duces to estimating the effective eigenvalue of the Markov 
matrix in terms of the magnetization autocorrelation 
function g(t) via A^*"* = g{t + l)/g{t). To estimate 
g{t) one would average over time products of the form 
77i(si)m(si+i). Eq. ^) would then yield g{t + 1) by 
replacing m{st) by the conditional expectation value of 
the magnetization at time t + 1, evaluated explicitly as 
Es^+i m{st+i)P{st+i\st). 

The crux is that the estimator of A^ satisfies a zero- 
variance principle, since Eqs. and contain an op- 
timizable trial state |/). In the ideal case, |/) is an exact 
eigenstate of the symmetrized Markov matrix P, and the 
'configurational eigenvalue' Al(s) equals the eigenvalue 
independent of s. Then, the estimator of the effective 
eigenvalue A^'' yields the exact eigenvalue without sta- 
tistical and systematic errors at finite t, if care is taken 
to arrange cancellation of the fluctuating factors in the 



estimators of H^*^ and N^*\ In practice, |/) is not an 
exact eigenstate, and this introduces statistical and sys- 
tematic errors. However, these errors are kept small by 
the zero- variance principle, if the trial states are accurate. 

Such optimized trial states are constructed prior to the 
main Monte Carlo run, by minimization of the variance 
of the configurational eigenvalue 

x'ip) = {ip-{p)fr)f. (9) 

As indicated, the variance depends on the parameters p 
of the trial state. Optimization over p is done following 
Umrigar et al. ]lO| : one samples M configurations Si, 
typically a few thousand, with probability ip'^Z~^ and 
approximates (p) by 

2. ^ EfU[fi^^,P)/Ms^)?[XL(.S.,p) - ^Ljp)? 

(10) 

Here Al denotes the weighted average of the configura- 
tional eigenvalue over the sample, while the modified no- 
tation explicitly shows dependences on the parameters p 
of the trial state |/). Near-optimal values of the parame- 
ters p can be obtained by minimization of the expression 
on the right-hand side of Eq. (|l^) for a fixed sample. 

Analysis of the exact left eigenvectors of the Markov 
matrix P for systems with L < 5 shows that the el- 
ements depend only on the magnetization to good ap- 
proximation. This suggests trial functions depending on 
long-wavelength components of the Fourier transform of 
Si, the zero-momentum component of which is just the 
magnetization m. The form 

/(s)=V^b(s) v^-^s), (11) 

where -^f^) — » ±ip'--^'> under spin inversion, yields an an- 
tisymmetric trial function, as required. The tilde in tps 
indicates that the temperature is used as a variational 
parameter, but we found that its optimal value is vir- 
tually indistinguishable from the true temperature. The 
^(±) 

were chosen as 
V'^^^ = E ak(m^)TOj^+^ -I- m E &k(?n^)™lr^ (12) 

k k 

^(-) = mEck(TO^)™L^' +E^k(TO^)m|^"\ (13) 

k k 

where the index k runs through a small set of multiplets 
of four or less long-wavelength wave vectors defining the 
m^^^ , translation and rotation symmetric sums of prod- 
ucts of Fourier transforms of the local magnetization; the 
k are selected so that is odd and is even un- 

der spin inversion; the coefficients ak, 6k, Ck and dk are 
polynomials of second order or less in m^. The degrees 
of these polynomials were chosen so that no terms oc- 
cur of higher degree than four in the spin variables. No 
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more than forty parameters were optimized for the trial 
functions used in the computations reported here. 

Since the probabihty distribution Eq. (^) is precisely 
the one purportedly generated by standard Monte Carlo, 
the sampling procedure is straightforward. The Monte 
Carlo algorithm used a random-number generator of the 
shift-register type. It was selected with care to avoid the 
introduction of systematic errors; see discussion and ref- 
erences m Ref. We used two Kirkpatrick-StoU 
generators, the results of which were combined by a bit- 
wise exclusive-or . For test purposes we replaced one 
Kirkpatrick-StoU generator by a linear congruential rule, 
but this did not reveal clear differences . 

For each system size 4 < L < 15, Monte Carlo averages 
were taken over 8 x 10* spin configurations. For L = 13 - 
15 these were separated by intervals of 16 sweeps (Monte 
Carlo steps per spin); 8 sweeps for L = 11 and 12; 2 
sweeps for L ~ 5 and 6; and only one sweep for L = A. 
The simulations of the remaining system sizes consisted 
of parts using intervals of 2, 4 or 8 sweeps. 

The numerical results for the effective second largest 
eigenvalue A^'' as a function of the projection time t ap- 
peared to converge rapidly. In agreement with scaled 
results for L < 5 spectra, we observe that convergence 
occurs within a few intervals as given above. Monte 
Carlo estimates of Al are shown in Table |, as are ex- 
act results for small systems. For system sizes L = 4 
and 5, the two types of calculation agree satisfactorily. 
The small numerical errors indicate that the variance- 
reducing method introduced above is quite effective. 

For finite system size L there are corrections to the 
leading scaling behavior tl . In the two-dimensional 
Ising model corrections to static equilibrium quantities 
occur with even powers of \/L |14|; thus we expect 
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where the series was arbitrarily truncated at order nc, 
but other powers of 1/L might occur as well. Ignoring 
the latter, we fitted the autocorrelation times of Table || 
to this form. Typical results of such fits are given in Ta- 
ble ||. The smallest systems do not fit Eq. ( 14 ) well, at 
least not for the ric values used. The residuals decrease 
rapidly when the minimum system size is increased and 
the consistency between the results for different sug- 
gests that Eq. (|lj) captures the essential scaling behavior 
of tl. From these results we chose the entry for L > 5 
and Tic = 2 as our best estimate: z = 2.1665 ± 0.0012, 
where we conservatively quote a two-sigma error. To our 
knowledge, this is the most precise estimate of z obtained 
to date, as evidenced by recent results summarized in Ta- 
ble III , The table shows that that the mutual consistency 
of the results for z has tended to improve in recent years. 
The only recent result that appears inconsistent with ours 
is due to Li et al. p^] . Its error is copied from Table I of 



Li et al.. The data in this table display finite-size depen- 
dences that exceed the quoted errors, which may explain 
the discrepancy with our result. 
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TABLE I. Second-largest eigenvalue of the Markov ma- 
trix. The first column indicates the method: exact numerical 
or Monte Carlo. 



method 


L 


Ai 


error 


exact 


2 


0.985702260395516 


0.000000000001 


exact 


3 


0.997409385126011 


0.000000000001 


exact 


4 


0.999245567376453 


0.000000000001 


exact 


5 


0.999708953624452 


0.000000000001 


MC 


4 


0.9992455685 


0.0000000094 


MC 


5 


0.9997089453 


0.0000000060 


MC 


6 


0.9998657194 


0.0000000045 


MC 


7 


0.9999299708 


0.0000000031 


MC 


8 


0.9999600854 


0.0000000023 


MC 


9 


0.9999756630 


0.0000000017 


MC 


10 


0.9999843577 


0.0000000014 


MC 


11 


0.9999895056 


0.0000000010 


MC 


12 


0.9999927107 


0.0000000008 


MC 


13 


0.9999947840 


0.0000000006 


MC 


14 


0.9999961736 


0.0000000005 


MC 


15 


0.9999971314 


0.0000000005 



TABLE II. Results of least-squares fits for the dynamic 
exponent. The first column shows the minimum system size 
included, the second the number of correction terms included, 
and the third column whether (y) or not (n) numerical exact 
results (for L < 5) are included. The last column contains 
the chi-square confidence index. ]l6|. 



L > 


nc 


exact 


z 


error 


Q 


4 


1 


n 


2.1769 


0.0001 


0.00 


5 


1 


n 


2.1705 


0.0002 


0.00 


6 


1 


n 


2.1688 


0.0003 


0.23 


7 


I 




9 1 fi7Q 


n onnfi 

u.uuuu 


n A-'\ 


8 


1 


n 


2.1672 


0.0010 


0.42 


4 


2 


n 


2.1650 


0.0003 


0.17 


5 


2 


n 


2.1665 


0.0006 


0.70 


6 


2 


n 


2.1662 


0.0013 


0.60 


7 


2 


n 


2.1648 


0.0024 


0.52 


4 


3 


n 


2.1672 


0.0009 


0.64 


5 


3 


n 


2.1656 


0.0020 


0.61 


6 


3 


n 


2.1625 


0.0044 


0.56 


3 


3 


y 


2.1653 


0.0004 


0.34 


4 


3 


y 


2.1670 


0.0009 


0.64 


5 


3 


y 


2.1657 


0.0020 


0.49 



TABLE HI. Comparison of recent results for the dynamic 
exponent z. Numerical errors are in parentheses. 



Reference 


Year 


Value 


Present work 


1996 


2.1665 (12) 


Li et al. |l| 


1995 


2.1337 (41) 


Linke et al. |l7|] 


1995 


2.160 (5) 


Grassberger [p^[ 


1995 


2.172 (6) 


Wang et al. ]19| 


1995 


2.16 (4) 


Baker and Erpenback [^ 


1994 


2.17 (1) 


Ito @ 


1993 


2.165 (10) 


Dammann and Reger [E2| 


1993 


2.183 (5) 


Matz et al. [||l 


1993 


2.35 (5) 


Miinkel et aL|^ 


1993 


2.21 (3) 


Stauffer || 


1993 


2.06 (2) 
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